Времянные ряды
неделя 4
Подключаем диск
In [2]:
from google.colab import drive
drive.mount('/content/gdrive',force_remount=True)
Переходим в директорию, загружаем библиотеки
In [3]:
cd /content/gdrive/'My Drive'/Colab_Notebooks/ML_c6_TAXI
In [0]:
PATH_TO_DATA = 'data'
In [0]:
import warnings
warnings.filterwarnings('ignore')
from sklearn.base import BaseEstimator, TransformerMixin
import os
import pandas as pd
import numpy as np
from scipy import stats
from tqdm import tqdm_notebook
from matplotlib import pyplot as plt
import seaborn as sns
sns.set_style("darkgrid")
import gc
import statsmodels.api as sm
from sklearn.linear_model import LinearRegression
from itertools import product
from sklearn.metrics import mean_absolute_error as mae
from sklearn.cluster import KMeans, AgglomerativeClustering
import scipy.cluster.hierarchy as shc
from scipy.cluster.hierarchy import dendrogram
from yellowbrick.cluster import KElbowVisualizer, SilhouetteVisualizer
from sklearn.manifold import TSNE,MDS
import folium
from sklearn.preprocessing import PowerTransformer
from scipy import stats
from dateutil.relativedelta import relativedelta
Список рассматриваемых зон
In [6]:
zones = ['1075', '1076', '1077', '1125', '1126', '1127', '1128', '1129', '1130', '1131', '1132', '1172', '1173', '1174', '1175', \
'1176', '1177', '1178', '1179', '1180', '1181', '1182', '1183', '1184', '1221', '1222', '1223', '1224', '1225', '1227',\
'1228', '1229', '1230', '1231', '1232', '1233', '1234', '1235', '1272', '1273', '1274', '1278', '1279', '1280', '1281',\
'1282', '1283', '1284', '1285', '1286', '1287', '1326', '1327', '1331', '1332', '1333', '1334', '1335', '1336', '1337', \
'1338', '1339', '1376', '1377', '1378', '1380', '1382', '1383', '1384', '1385', '1386', '1387', '1388', '1389', '1390', \
'1426', '1431', '1434', '1435', '1436', '1437', '1438', '1439', '1441', '1442', '1480', '1482', '1483', '1530', '1532', \
'1533', '1580', '1630', '1684', '1733', '1734', '1783', '2068', '2069', '2118', '2119', '2168']
zones = list(map(np.int64,zones))
len(zones)
Out[6]:
Загружаем 6 месяцев (декабрь 2015 г. - март 2016 г.)
In [7]:
%%time
time_series_6_month = pd.read_csv(os.path.join(PATH_TO_DATA,'time_series_6_month.csv'),index_col=0)
time_series_6_month.columns = range(1,2501)
dt = [pd.to_datetime(tt) for tt in list(time_series_6_month.index)]
time_series_6_month.index = dt
time_series_6_month.head()
Выделяем апрель, на нем будет проводится кластеризация
In [0]:
data_april = time_series_6_month.loc[pd.to_datetime('2016-04-01 00:00:00'):pd.to_datetime('2016-04-30 23:00:00')].copy()
In [9]:
data_april_zones = data_april[zones].copy()
data_april_zones.head()
Out[9]:
Визуальный анализ рядов
In [10]:
def time_series_plot(ts,zone_id):
with plt.xkcd():
plt.figure(figsize=(50,6))
plt.plot(ts, color = 'blue',linewidth = 1.5)
plt.grid(True)
plt.xlabel('time')
plt.ylabel('trip count')
plt.title(f'Zone {zone_id}')
for col in data_april_zones.columns[:2]:
time_series_plot(data_april_zones[col],col)
Нормализация
In [0]:
def normalisation(df):
normed_df = df.copy()
norm_coef = {}
for col in df.columns:
norm_coef[col]=[normed_df[col].mean(),normed_df[col].std()]
normed_df[col] = (normed_df[col]-norm_coef[col][0])/norm_coef[col][1]
return normed_df,norm_coef
In [12]:
data_zones_april_normed,norm_coef = normalisation(data_april_zones)
for col in data_zones_april_normed.columns[:2]:
time_series_plot(data_zones_april_normed[col],col)
Создадим признаки для рядов на основе статистик и данных по автокорреляции
In [0]:
def create_fe(data):
df = data.copy()
features={}
for col in data.columns:
#автокорреляционные признаки
fe_list0 = list(sm.graphics.tsa.acf(df[col].values.squeeze(),nlags=168))
fe_list = fe_list0[1:]
#медиана
fe_list.append(df[col].median())
#признаки, основанные на количестве пиков и впадин
row_shifted = df[col] - df[col].min()
fe_list.append(sum(row_shifted>0.9*row_shifted.max())*720/row_shifted.shape[0])
fe_list.append(sum(row_shifted<0.1*row_shifted.max())*720/row_shifted.shape[0])
fe_list.append((df[col].max()-df[col].min())/df[col].max())
#индикатор основной активности (выходные или рабочие дни)
st_we = [st for st in list(df[col].index) if st.weekday() == 5 or st.weekday() == 6]
we_count = sum(df[col].loc[st_we])
st_wd = [st for st in list(df[col].index) if st.weekday() != 5 and st.weekday() != 6]
wd_count = sum(df[col].loc[st_wd])
w_w_indic = int(2.5*we_count / wd_count > 1)
fe_list.append(w_w_indic)
features[col] = fe_list
feats_df = pd.DataFrame(features).T
return feats_df
Попробуем воспользоваться kmeans и примерно подобрать число кластеров с учетом качества и скорости решения.
In [14]:
features = create_fe(data_zones_april_normed)
model01 = KMeans(random_state=42)
visualizer = KElbowVisualizer(model01, k=(4,20), timings=False)
visualizer.fit(features)
visualizer.poof()
In [15]:
model02 = KMeans(10,random_state=42)
visualizer = SilhouetteVisualizer(model02)
visualizer.fit(features)
visualizer.poof()
Предварительно, 10 кластеров оптимально (см. график для правила локтя и коэффициента силуэта). Посмотрим на результаты кластеризации на признаках.
In [0]:
model1 = KMeans(n_clusters=10,random_state=42).fit(features)
In [17]:
def cluster_center(data):
data_to_work = data.copy().T
result = [data_to_work[column].mean() for column in data_to_work.columns]
return pd.Series(result,index=data_to_work.columns)
def plot_clusters(model_fitted,data):
def cluster_center(data):
data_to_work = data.copy().T
result = [data_to_work[column].mean() for column in data_to_work.columns]
return pd.Series(result,index=data_to_work.columns)
color_dict = {0:'red',1:'blue',2:'green',3:'purple',4:'orange',\
5:'darkred',6:'cadetblue',7:'darkblue',8:'black',\
9:'gray'}
frame_to_plot = data.T.copy()
frame_to_plot['cluster'] = model_fitted.labels_
overall_dist_list = []
for (c,sub_df) in frame_to_plot.groupby(by='cluster'):
plt.figure(figsize=(20,4))
cl_center = cluster_center(sub_df.T)
dist_list = []
for idx in sub_df.index:
dist = np.sqrt(sum([(coord - coord_c)**2 for (coord,coord_c) in zip(sub_df.loc[idx],cl_center)]))
dist_list.append(dist)
plt.plot(sub_df.loc[idx].values[:168],color = color_dict[c],alpha = 0.2)
mean_dist_c = np.mean(dist_list)
overall_dist_c = mean_dist_c*sub_df.shape[0]
overall_dist_list.append(overall_dist_c)
plt.plot(cl_center.values[:168],color = color_dict[c],linewidth=2)
plt.title(f'Кластер {c},число объектов в кластере {sub_df.shape[0]}, среднее отклонение от центра кластера {round(mean_dist_c,2)}');
print(f'суммарное отклонение от центров кластеров {round(sum(overall_dist_list),2)}')
print(f'среднее отклонение от центров кластеров {round(sum(overall_dist_list)/102,2)}')
plot_clusters(model1,data_zones_april_normed)
Видно, что некоторые кластеры имеют очень плохое качество, таким образом, информации в признаках мало. Попробуем воспользоваться агломеративной иерархической кластеризацией на значениях рядов за апрель.
In [0]:
model2 = AgglomerativeClustering(n_clusters=10).fit(data_zones_april_normed.T)
In [19]:
plot_clusters(model2,data_zones_april_normed)
Видно, что эти 10 кластеров имеют гораздо лучшее качество. Посмотрим на дендрограмму и определим, нельзя ли сократить количество кластеров, поскольку подбор 10 моделей сравнительно сложен.
In [20]:
with plt.xkcd():
plt.figure(figsize=(16, 10))
plt.title("Taxi counts dendrogram", fontsize=22)
dend = dendrogram(shc.linkage(data_zones_april_normed.T, method='ward'), labels=zones, color_threshold=37)
plt.xticks(fontsize=10)
plt.show()
In [21]:
model3 = AgglomerativeClustering(n_clusters=1).fit(data_zones_april_normed.T)
plot_clusters(model3,data_zones_april_normed)
Видно, какое низкое качество дает 1 кластер.
In [22]:
model4 = AgglomerativeClustering(n_clusters=5).fit(data_zones_april_normed.T)
plot_clusters(model4,data_zones_april_normed)
Итак, 5 кластеров дают достаточное качество для первого приближения.
In [23]:
model_final = model4
model_final.labels_
Out[23]:
In [26]:
id_cluster_df = pd.DataFrame([zones,model_final.labels_]).T
id_cluster_df.head()
Out[26]:
In [24]:
regions = pd.read_csv(os.path.join(PATH_TO_DATA,'regions.csv'),sep=';')
regions.head()
Out[24]:
Построим кластеры на карте.
In [27]:
m = folium.Map(location=[40.748817,-73.985428,], zoom_start=10.5)
color_dict = {0:'red',1:'blue',2:'green',3:'purple',4:'orange',\
5:'darkred',6:'cadetblue',7:'darkblue',8:'black',\
9:'gray'}
for cluster_label in np.unique(model_final.labels_):
cluster_id_list = list(id_cluster_df[id_cluster_df[1] == cluster_label][0].values)
cluster_indices = [id - 1 for id in cluster_id_list]
cluster_regions = regions.iloc[cluster_indices].copy()
#Центры ячеек для размещения маркеров
cluster_regions['mid_long']=cluster_regions.apply(lambda row: (row['east']+row['west'])/2, axis=1)
cluster_regions['mid_lat']=cluster_regions.apply(lambda row: (row['north']+row['south'])/2, axis=1)
taxi_markers = folium.map.FeatureGroup(name=f'Cluster {cluster_label}')
for row in cluster_regions.values:
id = int(row[0])
message = f'<i>Id: {id}</i>'
coords = [row[-1],row[-2]]
taxi_markers.add_children(folium.Marker(coords,
popup=message,
icon=folium.Icon(icon='taxi',color=color_dict[cluster_label],prefix='fa')
)
)
taxi_markers.add_to(m)
folium.LayerControl().add_to(m)
m
Out[27]:
Видно, что логика в разбиении присутствует, можно заметить, что аэропорты объединились в 1 кластер, также можно заметить, что центр города и окраина принадлежат разным кластерам.
Построим характерные ряды для каждого кластера (нормализованные центры), подберем модели на данных до апреля, валдация будет на апреле. Оптимальные параметры запишем в словарь.
In [28]:
time_series_4_month = time_series_6_month.loc[:pd.to_datetime('2016-04-30 23:00:00')]
time_series_4_month.tail()
Out[28]:
In [29]:
ts_normed,norm_coef = normalisation(time_series_4_month[zones])
cluster_centers = {}
for cluster_label in np.unique(model_final.labels_):
cluster_id_list = list(id_cluster_df[id_cluster_df[1] == cluster_label][0].values)
cluster_centers[cluster_label] = cluster_center(ts_normed[cluster_id_list])
time_series_centers = pd.DataFrame(cluster_centers)
time_series_centers.tail()
Out[29]:
In [30]:
for col in time_series_centers.columns:
time_series_plot(time_series_centers[col].iloc[-720:].iloc[:168],col)
In [0]:
def reg_feats(start,stop,kw=2,ky=0,linear_fe=False):
#Вспомогательная функция для построения признаков
arg = np.array([x for x in range(start,stop+1)])
fe_list = []
columns = []
#Линейный признак для описания простейшего тренда
if linear_fe:
fe_list.append(arg)
columns.append('t')
#Признаки для описания недельной сезонности
for i in range(1,kw+1):
fe_list.append(np.sin(arg*2*np.pi*i/168))
columns.append(f'sin(2*pi*t*{i}/168)')
fe_list.append(np.cos(arg*2*np.pi*i/168))
columns.append(f'cos(2*pi*t*{i}/168)')
#Признаки для описания годовой сезонности
for j in range(1,ky+1):
fe_list.append(np.sin(arg*2*np.pi*j/8766))
columns.append(f'sin(2*pi*t*{j}/8766)')
fe_list.append(np.cos(arg*2*np.pi*i/8766))
columns.append(f'cos(2*pi*t*{j}/8766)')
df = pd.DataFrame(fe_list).T
df.columns = columns
return df
In [0]:
ts0 = time_series_centers[0]
ts = ts0.copy()
ts_work = ts
ts_tr = ts_work.iloc[:-720]
ts_val = ts_work.iloc[-720:]
feats = reg_feats(1,ts_tr.shape[0],kw=3,ky=0)
feats.head()
Out[0]:
In [0]:
def reg_res_get_plot(ts,reg_features,regr=LinearRegression()):
feats = reg_features
regressor = regr
regressor.fit(feats,ts)
preds = regressor.predict(feats)
plt.figure(figsize=(50,16))
ax=plt.subplot(211)
ax.plot(ts.index[-720:],preds[-720:])
plt.title('regression');
resids = ts - preds
ax=plt.subplot(212)
ax.plot(resids[-720:])
plt.title('regression residuals');
return resids
resids = reg_res_get_plot(ts_tr,reg_features=reg_feats(1,ts_tr.shape[0],kw=2,ky=0),regr=LinearRegression())
In [0]:
ts_df = resids.to_frame(name="resids")
def stl_plot(ts,lim):
with plt.xkcd():
result = sm.tsa.seasonal_decompose(ts)
plt.figure(figsize=(50,6))
result.observed[-lim:].plot(title = 'observed')
plt.grid(True)
plt.figure(figsize=(50,6))
result.trend[-lim:].plot(title = 'trend')
plt.grid(True)
plt.figure(figsize=(50,6))
result.seasonal[-lim:].plot(title = 'seasonal')
plt.grid(True)
plt.figure(figsize=(50,6))
result.resid[-lim:].plot(title = 'residuals')
plt.grid(True)
stl_plot(ts_df.resids,lim=500)
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(resids)[1])
In [0]:
ts_df['diff1_resids'] = ts_df.resids - ts_df.resids.shift(24)
stl_plot(ts_df.diff1_resids[25:],lim=500)
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(ts_df.diff1_resids[25:] )[1])
In [0]:
def acf_pacf_plot(series,n_lags):
plt.figure(figsize=(15,10))
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(series.values.squeeze(), lags=n_lags, ax=ax);
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(series.values.squeeze(), lags=n_lags, ax=ax);
acf_pacf_plot(ts_df.diff1_resids[25:],n_lags=168)
In [0]:
ps = range(0, 4)
d=0
qs = range(0, 4)
Ps = range(0, 2)
D=1
Qs = range(0, 3)
In [0]:
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)
Out[0]:
In [0]:
class TimeSeriesModel(BaseEstimator, TransformerMixin):
'''
Модель для предсказания временного ряда.
Использует регрессию на тригонометрические признаки и модель SARIMA
'''
def __init__(self,kw=2,ky=0,linear_fe=False,params=(1,1,1),seasonal_params=(1,1,1,24),use_yj=False):
self.kw = kw
self.ky = ky
self.linear_fe = linear_fe
self.params=params
self.seasonal_params=seasonal_params
self.use_yj = use_yj
def reg_feats(self,start,stop,kw=2,ky=0,linear_fe=False):
#Вспомогательная функция для построения признаков
arg = np.array([x for x in range(start,stop+1)])
fe_list = []
columns = []
#Линейный признак для описания простейшего тренда
if linear_fe:
fe_list.append(arg)
columns.append('t')
#Признаки для описания недельной сезонности
for i in range(1,kw+1):
fe_list.append(np.sin(arg*2*np.pi*i/168))
columns.append(f'sin(2*pi*t*{i}/168)')
fe_list.append(np.cos(arg*2*np.pi*i/168))
columns.append(f'cos(2*pi*t*{i}/168)')
#Признаки для описания годовой сезонности
for j in range(1,ky+1):
fe_list.append(np.sin(arg*2*np.pi*j/8766))
columns.append(f'sin(2*pi*t*{j}/8766)')
fe_list.append(np.cos(arg*2*np.pi*i/8766))
columns.append(f'cos(2*pi*t*{j}/8766)')
df = pd.DataFrame(fe_list).T
df.columns = columns
return df
def fit(self,ts,save_ar_params=False,load_ar_params=False,PATH_TO_MODELS=''):
#Обучение, вход - df
#Нормализация для сглаживания пиков и стабилизации дисперсии исходного ряда, если нужно
if self.use_yj:
self.yj = PowerTransformer(method='yeo-johnson', standardize=False, copy=True)
self.yj.fit(ts)
tr_ts = pd.DataFrame(self.yj.transform(ts.copy()), index=ts.index,columns=[ts.columns[0]])
self.ts = tr_ts
else:
self.ts = ts.copy()
#Признаки для регрессии и обучение регрессионной модели
feats = self.reg_feats(start=1,stop=self.ts.shape[0],kw=self.kw, ky=self.ky,linear_fe=self.linear_fe)
self.regressor = LinearRegression()
self.regressor.fit(feats,self.ts)
#Получаем остатки регрессии для использования в SARIMA
preds = self.regressor.predict(feats)
resids = self.ts - preds
#Загрузка или обучение авторегрессионной модели
if load_ar_params:
#Загрузка параметров
ar_params = pd.read_csv(os.path.join(PATH_TO_MODELS,f'{self.ts.columns[0]}_params'),\
index_col=0,names=['params'],header=0).params
ar_params.name=None
#фильтр (быстро)
self.model = sm.tsa.statespace.SARIMAX(resids, order=self.params,
seasonal_order=self.seasonal_params).filter(ar_params)
else:
#обучение (долго)
self.model=sm.tsa.statespace.SARIMAX(resids, order=self.params,
seasonal_order=self.seasonal_params).fit(disp=-1)
#Сохранение авторегрессионной модели
if not load_ar_params and save_ar_params:
self.model.params.to_csv(os.path.join(PATH_TO_MODELS,f'{self.ts.columns[0]}_params'))
return self
def extend(self,bigger_ts,refit_regressor = False):
#Расширяет модель на более широкий диапазон данных, не переобучая SARIMAX, вход - df
#Нормализация для сглаживания пиков и стабилизации дисперсии исходного ряда, если нужно
if self.use_yj:
self.yj = PowerTransformer(method='yeo-johnson', standardize=False, copy=True)
self.yj.fit(bigger_ts)
tr_ts = pd.DataFrame(self.yj.transform(bigger_ts.copy()), index=bigger_ts.index,columns=[bigger_ts.columns[0]])
self.ts = tr_ts
else:
self.ts = bigger_ts.copy()
feats = self.reg_feats(start=1,stop=self.ts.shape[0],kw=self.kw, ky=self.ky,linear_fe=self.linear_fe)
if refit_regressor:
self.regressor.fit(feats,self.ts)
preds = self.regressor.predict(feats)
resids = self.ts - preds
new_model = sm.tsa.statespace.SARIMAX(resids, order=self.params,
seasonal_order=self.seasonal_params).filter(self.model.params)
self.model = new_model
return self
def predict(self,start,stop,verbose=False):
#Предсказание, выход df
#start, stop - timestamps
#Получаем признаки регрессии для предсказания
#Timestamps to int
ts_index = list(self.ts.index)
if start in ts_index:
start_pos = ts_index.index(start) + 1
else:
start_pos = self.ts.shape[0] + (start - ts_index[len(ts_index)-1]).total_seconds()/3600 + 1
stop_pos = start_pos + (stop - start).total_seconds()/3600
feats = self.reg_feats(int(start_pos),int(stop_pos),self.kw,self.ky,self.linear_fe)
#Предсказание регрессии
pred_regression = self.regressor.predict(feats)
#Предсказание SARIMA
pred_resids = self.model.predict(start=start, end=stop)
#Итоговое предсказание как суперпозиция
result = pd.DataFrame(pred_resids) + pred_regression
result.columns = [f'{self.ts.columns[0]}_forecast']
if self.use_yj:
index = result.index
result = pd.DataFrame(self.yj.inverse_transform(result),index = index,columns=[f'{self.ts.columns[0]}_forecast'])
if verbose:
plt.figure()
plt.plot(pd.DataFrame(pred_resids,index = result.index),color='g')#график предсказанных остатков для отладки
plt.plot(pd.DataFrame(pred_regression,index = result.index),color='r') #график регрессии для отладки
return result
In [0]:
lim = 168
results = []
best_aic = float("inf")
warnings.filterwarnings('ignore')
for param in tqdm_notebook(parameters_list):
#try except нужен, потому что на некоторых наборах параметров модель не обучается
try:
model = sm.tsa.statespace.SARIMAX(resids[-lim:], order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], 24)).fit(disp=-1)
#выводим параметры, на которых модель не обучается и переходим к следующему набору
except ValueError:
print('wrong parameters:', param)
continue
aic = model.aic
#сохраняем лучшую модель, aic, параметры
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])
warnings.filterwarnings('default')
In [0]:
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True)[:10])
In [0]:
model = TimeSeriesModel(kw=2,ky=0,linear_fe=False,params=(best_param[0],d,best_param[1]),seasonal_params=(best_param[2],D,best_param[3],24),use_yj=False)
In [0]:
%%time
model.fit(pd.DataFrame(ts_tr[-720:]))
print(model.model.summary())
In [0]:
plt.figure(figsize=(15,8))
ax = plt.subplot(211)
ax.plot(model.model.resid[-720:])
plt.ylabel('Residuals')
ax = plt.subplot(212)
sm.graphics.tsa.plot_acf(model.model.resid[25:].squeeze(), lags=168, ax=ax)
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(model.model.resid[25:], 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(model.model.resid[25:])[1])
In [0]:
index=[pd.to_datetime(tt) for tt in ts_tr.index[-500:]]
ts_resids= pd.Series(model.model.resid[-500:],index=index)
stl_plot(ts_resids,lim=500)
In [0]:
prediction = model.predict(pd.to_datetime('2016-04-01 00:00:00'),pd.to_datetime('2016-04-30 23:00:00'),verbose=True)
In [0]:
plt.figure(figsize=(20,10))
plt.plot(ts_val,color='b')
plt.plot(ts_val.index,prediction,color='r');
In [0]:
error = mae(ts_val,prediction)
error
Out[0]:
In [0]:
'''
cluster0_params:
kw=2
SARIMAX(2, 0, 1)x(1, 1, 1, 24)
use_yj=False
'''
Out[0]:
In [0]:
cl0_model_params= {'kw':2,
'params':(2, 0, 1),
'seasonal_params':(1, 1, 1, 24),
'use_yj':False
}
In [0]:
ts1 = time_series_centers[1]
ts = ts1.copy()
yj = PowerTransformer(method='yeo-johnson', standardize=False, copy=True)
yj.fit(ts.values.reshape(-1,1))
ts_work = pd.Series(yj.transform(ts.values.reshape(-1,1)).reshape(-1),index=ts.index)
ts_tr = ts_work.iloc[:-720]
ts_val = ts_work.iloc[-720:]
feats = reg_feats(1,ts_tr.shape[0],kw=3,ky=0)
feats.head()
Out[0]:
In [0]:
resids = reg_res_get_plot(ts_tr,reg_features=reg_feats(1,ts_tr.shape[0],kw=3,ky=0),regr=LinearRegression())
In [0]:
ts_df = resids.to_frame(name="resids")
stl_plot(ts_df.resids,lim=500)
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(resids)[1])
In [0]:
ts_df['diff1_resids'] = ts_df.resids - ts_df.resids.shift(24)
stl_plot(ts_df.diff1_resids[25:],lim=500)
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(ts_df.diff1_resids[25:] )[1])
In [0]:
acf_pacf_plot(ts_df.diff1_resids[25:],n_lags=168)
In [0]:
ps = range(0, 2)
d=0
qs = range(0, 4)
Ps = range(0, 4)
D=1
Qs = range(0, 2)
In [0]:
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)
Out[0]:
In [0]:
lim = 168
results = []
best_aic = float("inf")
warnings.filterwarnings('ignore')
for param in tqdm_notebook(parameters_list):
#try except нужен, потому что на некоторых наборах параметров модель не обучается
try:
model = sm.tsa.statespace.SARIMAX(resids[-lim:], order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], 24)).fit(disp=-1)
#выводим параметры, на которых модель не обучается и переходим к следующему набору
except ValueError:
print('wrong parameters:', param)
continue
aic = model.aic
#сохраняем лучшую модель, aic, параметры
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])
warnings.filterwarnings('default')
In [0]:
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True)[:10])
In [0]:
model = TimeSeriesModel(kw=3,ky=0,linear_fe=False,params=(best_param[0],d,best_param[1]),seasonal_params=(best_param[2],D,best_param[3],24),use_yj=True)
In [0]:
%%time
model.fit(pd.DataFrame(ts_tr[-720:]))
print(model.model.summary())
In [0]:
plt.figure(figsize=(15,8))
ax = plt.subplot(211)
ax.plot(model.model.resid[-720:])
plt.ylabel('Residuals')
ax = plt.subplot(212)
sm.graphics.tsa.plot_acf(model.model.resid[25:].squeeze(), lags=168, ax=ax)
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(model.model.resid[25:], 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(model.model.resid[25:])[1])
In [0]:
index=[pd.to_datetime(tt) for tt in ts_tr.index[-500:]]
ts_resids= pd.Series(model.model.resid[-500:],index=index)
stl_plot(ts_resids,lim=500)
In [0]:
prediction = model.predict(pd.to_datetime('2016-04-01 00:00:00'),pd.to_datetime('2016-04-30 23:00:00'),verbose=True)
In [0]:
plt.figure(figsize=(20,10))
plt.plot(ts_val,color='b')
plt.plot(ts_val.index,prediction,color='r');
In [0]:
error = mae(ts_val,prediction)
error
Out[0]:
In [0]:
'''
cluster1_params:
kw=3
SARIMAX(0, 0, 3)x(0, 1, 1, 24)
use_yj=True
'''
Out[0]:
In [0]:
cl1_model_params= {'kw':3,
'params':(0, 0, 3),
'seasonal_params':(0, 1, 1, 24),
'use_yj':True
}
In [0]:
ts2 = time_series_centers[2]
ts = ts2.copy()
yj = PowerTransformer(method='yeo-johnson', standardize=False, copy=True)
yj.fit(ts.values.reshape(-1,1))
ts_work = pd.Series(yj.transform(ts.values.reshape(-1,1)).reshape(-1),index=ts.index)
ts_tr = ts_work.iloc[:-720]
ts_val = ts_work.iloc[-720:]
feats = reg_feats(1,ts_tr.shape[0],kw=3,ky=0)
feats.head()
Out[0]:
In [0]:
resids = reg_res_get_plot(ts_tr,reg_features=reg_feats(1,ts_tr.shape[0],kw=2,ky=0),regr=LinearRegression())
In [0]:
ts_df = resids.to_frame(name="resids")
stl_plot(ts_df.resids,lim=500)
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(resids)[1])
In [0]:
ts_df['diff1_resids'] = ts_df.resids - ts_df.resids.shift(24)
stl_plot(ts_df.diff1_resids[25:],lim=500)
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(ts_df.diff1_resids[25:] )[1])
In [0]:
acf_pacf_plot(ts_df.diff1_resids[25:],n_lags=168)
In [0]:
ps = range(0, 3)
d=0
qs = range(0, 4)
Ps = range(0, 2)
D=1
Qs = range(0, 3)
In [0]:
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)
Out[0]:
In [0]:
lim = 168
results = []
best_aic = float("inf")
warnings.filterwarnings('ignore')
for param in tqdm_notebook(parameters_list):
#try except нужен, потому что на некоторых наборах параметров модель не обучается
try:
model = sm.tsa.statespace.SARIMAX(resids[-lim:], order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], 24)).fit(disp=-1)
#выводим параметры, на которых модель не обучается и переходим к следующему набору
except ValueError:
print('wrong parameters:', param)
continue
aic = model.aic
#сохраняем лучшую модель, aic, параметры
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])
warnings.filterwarnings('default')
In [0]:
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True)[:10])
In [0]:
model = TimeSeriesModel(kw=2,ky=0,linear_fe=False,params=(best_param[0],d,best_param[1]),seasonal_params=(best_param[2],D,best_param[3],24),use_yj=True)
In [0]:
%%time
model.fit(pd.DataFrame(ts_tr[-720:]))
print(model.model.summary())
In [0]:
plt.figure(figsize=(15,8))
ax = plt.subplot(211)
ax.plot(model.model.resid[-720:])
plt.ylabel('Residuals')
ax = plt.subplot(212)
sm.graphics.tsa.plot_acf(model.model.resid[25:].squeeze(), lags=168, ax=ax)
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(model.model.resid[25:], 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(model.model.resid[25:])[1])
In [0]:
index=[pd.to_datetime(tt) for tt in ts_tr.index[-500:]]
ts_resids= pd.Series(model.model.resid[-500:],index=index)
stl_plot(ts_resids,lim=500)
In [0]:
prediction = model.predict(pd.to_datetime('2016-04-01 00:00:00'),pd.to_datetime('2016-04-30 23:00:00'),verbose=True)
In [0]:
plt.figure(figsize=(20,10))
plt.plot(ts_val,color='b')
plt.plot(ts_val.index,prediction,color='r');
In [0]:
error = mae(ts_val,prediction)
error
Out[0]:
In [0]:
'''
cluster2_params:
kw=2
SARIMAX(2, 0, 1)x(1, 1, 1, 24)
use_yj=True
'''
Out[0]:
In [0]:
cl2_model_params= {'kw':2,
'params':(2, 0, 1),
'seasonal_params':(1, 1, 1, 24),
'use_yj':True
}
In [0]:
ts3 = time_series_centers[3]
ts = ts3.copy()
ts_work = ts
ts_tr = ts_work.iloc[:-720]
ts_val = ts_work.iloc[-720:]
feats = reg_feats(1,ts_tr.shape[0],kw=3,ky=0)
feats.head()
Out[0]:
In [0]:
resids = reg_res_get_plot(ts_tr,reg_features=reg_feats(1,ts_tr.shape[0],kw=2,ky=0),regr=LinearRegression())
In [0]:
ts_df = resids.to_frame(name="resids")
stl_plot(ts_df.resids,lim=500)
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(resids)[1])
In [0]:
ts_df['diff1_resids'] = ts_df.resids - ts_df.resids.shift(24)
stl_plot(ts_df.diff1_resids[25:],lim=500)
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(ts_df.diff1_resids[25:] )[1])
In [0]:
acf_pacf_plot(ts_df.diff1_resids[25:],n_lags=168)
In [0]:
ps = range(0, 3)
d=0
qs = range(0, 4)
Ps = range(0, 2)
D=1
Qs = range(0, 4)
In [0]:
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)
Out[0]:
In [0]:
lim = 168
results = []
best_aic = float("inf")
warnings.filterwarnings('ignore')
for param in tqdm_notebook(parameters_list):
#try except нужен, потому что на некоторых наборах параметров модель не обучается
try:
model = sm.tsa.statespace.SARIMAX(resids[-lim:], order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], 24)).fit(disp=-1)
#выводим параметры, на которых модель не обучается и переходим к следующему набору
except ValueError:
print('wrong parameters:', param)
continue
aic = model.aic
#сохраняем лучшую модель, aic, параметры
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])
warnings.filterwarnings('default')
In [0]:
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True)[:10])
In [0]:
model = TimeSeriesModel(kw=2,ky=0,linear_fe=False,params=(best_param[0],d,best_param[1]),seasonal_params=(best_param[2],D,best_param[3],24),use_yj=False)
In [0]:
%%time
model.fit(pd.DataFrame(ts_tr[-720:]))
print(model.model.summary())
In [0]:
plt.figure(figsize=(15,8))
ax = plt.subplot(211)
ax.plot(model.model.resid[-720:])
plt.ylabel('Residuals')
ax = plt.subplot(212)
sm.graphics.tsa.plot_acf(model.model.resid[25:].squeeze(), lags=168, ax=ax)
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(model.model.resid[25:], 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(model.model.resid[25:])[1])
In [0]:
index=[pd.to_datetime(tt) for tt in ts_tr.index[-500:]]
ts_resids= pd.Series(model.model.resid[-500:],index=index)
stl_plot(ts_resids,lim=500)
In [0]:
prediction = model.predict(pd.to_datetime('2016-04-01 00:00:00'),pd.to_datetime('2016-04-30 23:00:00'),verbose=True)
In [0]:
plt.figure(figsize=(20,10))
plt.plot(ts_val,color='b')
plt.plot(ts_val.index,prediction,color='r');
In [0]:
error = mae(ts_val,prediction)
error
Out[0]:
In [0]:
'''
cluster3_params:
kw=2
SARIMAX(2, 0, 2)x(1, 1, 1, 24)
use_yj=False
'''
Out[0]:
In [0]:
cl3_model_params= {'kw':2,
'params':(2, 0, 2),
'seasonal_params':(1, 1, 1, 24),
'use_yj':False
}
In [0]:
ts4 = time_series_centers[4]
ts = ts4.copy()
yj = PowerTransformer(method='yeo-johnson', standardize=False, copy=True)
yj.fit(ts.values.reshape(-1,1))
ts_work = pd.Series(yj.transform(ts.values.reshape(-1,1)).reshape(-1),index=ts.index)
ts_tr = ts_work.iloc[:-720]
ts_val = ts_work.iloc[-720:]
feats = reg_feats(1,ts_tr.shape[0],kw=3,ky=0)
feats.head()
Out[0]:
In [0]:
resids = reg_res_get_plot(ts_tr,reg_features=reg_feats(1,ts_tr.shape[0],kw=3,ky=0),regr=LinearRegression())
In [0]:
ts_df = resids.to_frame(name="resids")
stl_plot(ts_df.resids,lim=500)
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(resids)[1])
In [0]:
ts_df['diff1_resids'] = ts_df.resids - ts_df.resids.shift(24)
stl_plot(ts_df.diff1_resids[25:],lim=500)
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(ts_df.diff1_resids[25:] )[1])
In [0]:
acf_pacf_plot(ts_df.diff1_resids[25:],n_lags=168)
In [0]:
ps = range(0, 3)
d=0
qs = range(0, 4)
Ps = range(0, 3)
D=1
Qs = range(0, 3)
In [0]:
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)
Out[0]:
In [0]:
lim = 168
results = []
best_aic = float("inf")
warnings.filterwarnings('ignore')
for param in tqdm_notebook(parameters_list):
#try except нужен, потому что на некоторых наборах параметров модель не обучается
try:
model = sm.tsa.statespace.SARIMAX(resids[-lim:], order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], 24)).fit(disp=-1)
#выводим параметры, на которых модель не обучается и переходим к следующему набору
except ValueError:
print('wrong parameters:', param)
continue
aic = model.aic
#сохраняем лучшую модель, aic, параметры
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])
warnings.filterwarnings('default')
In [0]:
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True)[:10])
In [0]:
model = TimeSeriesModel(kw=3,ky=0,linear_fe=False,params=(best_param[0],d,best_param[1]),seasonal_params=(best_param[2],D,best_param[3],24),use_yj=True)
In [0]:
%%time
model.fit(pd.DataFrame(ts_tr[-720:]))
print(model.model.summary())
In [0]:
plt.figure(figsize=(15,8))
ax = plt.subplot(211)
ax.plot(model.model.resid[-720:])
plt.ylabel('Residuals')
ax = plt.subplot(212)
sm.graphics.tsa.plot_acf(model.model.resid[25:].squeeze(), lags=168, ax=ax)
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(model.model.resid[25:], 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(model.model.resid[25:])[1])
In [0]:
index=[pd.to_datetime(tt) for tt in ts_tr.index[-500:]]
ts_resids= pd.Series(model.model.resid[-500:],index=index)
stl_plot(ts_resids,lim=500)
In [0]:
prediction = model.predict(pd.to_datetime('2016-04-01 00:00:00'),pd.to_datetime('2016-04-30 23:00:00'),verbose=True)
In [0]:
plt.figure(figsize=(20,10))
plt.plot(ts_val,color='b')
plt.plot(ts_val.index,prediction,color='r');
In [0]:
error = mae(ts_val,prediction)
error
Out[0]:
In [0]:
'''
cluster4_params:
kw=3
SARIMAX(2, 0, 1)x(1, 1, 1, 24)
use_yj=True
'''
Out[0]:
In [0]:
cl4_model_params= {'kw':3,
'params':(2, 0, 1),
'seasonal_params':(1, 1, 1, 24),
'use_yj':True
}
Словарь с параметрами моделей для кластеров
In [37]:
cluster_params = [cl0_model_params,cl1_model_params,cl2_model_params,cl3_model_params,cl4_model_params]
cluster_params
Out[37]:
In [38]:
id_cluster_df
Out[38]:
In [0]:
model_params_dict = {}
for zone_info in id_cluster_df.values:
zone_id = zone_info[0]
zone_cluster_id = zone_info[1]
model_params = cluster_params[zone_cluster_id]
model_params_dict[zone_id] = model_params
Обучение моделей, сохранение параметров в файл. Обучаем на данных за декабрь.
In [0]:
path_to_params = os.path.join(PATH_TO_DATA,'models_params')
In [42]:
%%time
warnings.filterwarnings("ignore")
for zone in model_params_dict:
ts = time_series_6_month[[zone]]
model_params = model_params_dict[zone]
model = TimeSeriesModel(kw=model_params['kw'],ky=0,linear_fe=False,params=model_params['params'],\
seasonal_params=model_params['seasonal_params'],use_yj=model_params['use_yj'])
model.fit(ts.iloc[:720],save_ar_params=True,PATH_TO_MODELS=path_to_params)
Расширяем модели на большие диапазоны, предсказаваем май и июнь.
In [44]:
%%time
warnings.filterwarnings("ignore")
predictions_may_dict = {}
predictions_june_dict = {}
for zone in model_params_dict:
ts = time_series_6_month[[zone]]
model_params = model_params_dict[zone]
model = TimeSeriesModel(kw=model_params['kw'],ky=0,linear_fe=False,params=model_params['params'],\
seasonal_params=model_params['seasonal_params'],use_yj=model_params['use_yj'])
model.fit(ts.iloc[:720],load_ar_params=True,PATH_TO_MODELS=path_to_params)
model.extend(ts.loc[:pd.to_datetime('2016-04-30 23:00:00')],refit_regressor=True)
prediction_may = model.predict(pd.to_datetime('2016-05-01 00:00:00'),pd.to_datetime('2016-05-31 23:00:00'))
predictions_may_dict[zone] = prediction_may
model.extend(ts,refit_regressor=True)
prediction_june = model.predict(pd.to_datetime('2016-06-01 00:00:00'),pd.to_datetime('2016-06-30 23:00:00'))
predictions_june_dict[zone] = prediction_june
Сохраняем предсказанные ряды.
In [58]:
for zone in predictions_june_dict:
pred_june_df = pd.concat([predictions_june_dict[zone] for zone in predictions_june_dict],axis=1)
pred_june_df.columns = predictions_june_dict.keys()
pred_june_df.head()
Out[58]:
In [59]:
for zone in predictions_may_dict:
pred_may_df = pd.concat([predictions_may_dict[zone] for zone in predictions_may_dict],axis=1)
pred_may_df.columns = predictions_may_dict.keys()
pred_may_df.head()
Out[59]:
In [0]:
pred_june_df.to_csv('pred_june.csv')
pred_may_df.to_csv('pred_may.csv')
In [70]:
df_may = pd.read_csv('pred_may.csv',index_col = 0)
cols_may = [int(col) for col in df_may.columns]
df_may.columns = cols_may
dt_may = [pd.to_datetime(tt) for tt in list(df_may.index)]
df_may.index = dt_may
df_may.tail()
Out[70]:
In [71]:
df_june = pd.read_csv('pred_june.csv',index_col = 0)
cols_june = [int(col) for col in df_june.columns]
df_june.columns = cols_june
dt_june = [pd.to_datetime(tt) for tt in list(df_june.index)]
df_june.index = dt_june
df_june.tail()
Out[71]:
Ниде приведен набросок онлайн алгоритма, который можно было бы использовать для предсказания на 6 часов вперед. Расширяем модель до временной точки, затем предсказываем на 6 часов, записываем ответ в результат.
In [0]:
'''
%%time
time_points = list(ts.loc[pd.to_datetime('2016-04-30 23:00:00'):pd.to_datetime('2016-05-31 17:00:00')].index)
result_df = pd.DataFrame(columns=['y'])
for zone in model_dict:
ts = time_series_6_month[[zone]]
for time in time_points:
model_dict[zone].extend(ts.loc[:time],refit_regressor=True)
start = time + relativedelta(hours=1)
stop = time + relativedelta(hours=6)
prediction = model_dict[zone].predict(start,stop)
time_string = '_'.join([str(zone),str(time.date()),str(time.hour)])
for ind,value in enumerate(prediction.values):
val_id = '_'.join([time_string,str(ind+1)])
result_df.loc[val_id] = value
result_df
'''
В нашем случае, в связи с ограничениями во времени и ресурсах, воспользуемся менее точным, но гораздо более быстрым алгоритмом. Воспользуемся предсказанными уже данными в предыдущих пунктах и распределим их по временным точкам.
In [0]:
def create_result(df,time_points):
time_inds = []
for zone in zones:
for time_point in time_points:
for count in range(1,7):
time_string = '_'.join([str(zone),str(time_point.date()),str(time_point.hour),str(count)])
time_inds.append(time_string)
parts_list = []
for time in time_points:
part = df.loc[time+relativedelta(hours=1):time+relativedelta(hours=6)].copy()
parts_list.append(part)
values = pd.concat(parts_list).values.T.reshape(-1,1)
result = pd.DataFrame(values,index = time_inds,columns=['y'])
return result
Результат для мая
In [109]:
result_may = create_result(df_may,time_points)
result_may.tail()
Out[109]:
истинные данные за май
In [111]:
real_may = create_result(time_series_6_month[zones],time_points)
real_may.head()
Out[111]:
Ошибка предсказания мая
In [112]:
error_pred = mae(real_may,result_may)
error_pred
Out[112]:
Получим предсказания для июня, проверим на Kaggle
In [0]:
time_points_june = list(pd.date_range(start=pd.to_datetime('2016-05-31 23:00:00'),end=pd.to_datetime('2016-06-30 17:00:00'),freq='1H'))
answer_df = create_result(df_june,time_points_june)
In [0]:
answer_df.to_csv('submission.csv',index_label='id')
Итак, для модели, которая работала всего 1,5 часа на среднем железе такой результат, на мой взгляд, довольно неплох. Для улучшения результата можно подбирать модели не для 5 кластеров, а для 10, при этом вместо того, чтобы обучать на 1 месяце и расширять на больший диапазон, можно обучать сразу на всех данных.
In [0]: